#include "Vertex.h"

Eigenstructure* Vertex::eigenstructures = 0;

void Vertex::useEigenstructures(Eigenstructure* e) {
    eigenstructures = e;
}

Vertex::Vertex() {
    index = -1;
    N = 0;
    sharpness = 0;
    ccVertex = 0;
    opt = false;
}

Vertex::Vertex(int index, const Vector3d& pos, int sharpness, bool opt) {
    this->index = index;
    this->sharpness = sharpness;
    this->pos = pos;
    this->opt = opt;
    N = 0;
    ccVertex = 0;
}

void Vertex::sort() {
    vector<Edge*> sortedEdges;
    vector<Face*> sortedFaces;

    Face* currFace = faces[0];
    Vertex* currVertex = currFace->getFromVertex(this, -1);
    Edge* currEdge = getEdge(this->index, currVertex->index);
    do {
        sortedEdges.push_back(currEdge);
        sortedFaces.push_back(currFace);
        currVertex = currFace->getFromVertex(currVertex, 2);
        currEdge = getEdge(this->index, currVertex->index);
        currFace = currEdge->getRightFace(this);
    } while (currFace != faces[0]);

    edges = sortedEdges;
    faces = sortedFaces;
}

void Vertex::setCCVertex(Vertex* v) {
    if (sharpness) {
        ccVertex = v;
        ccVertex->pos = pos;
    } else {
        int sharpValence = 0;
        for (int i = 0; i < N; i++) {
            if (edges[i]->sharpness) sharpValence++;
        }
        Vector3d nextPos;
        Vector3d faceAverage = Vector3d(0, 0, 0);
        Vector3d edgeAverage = Vector3d(0, 0, 0);
        
        switch (sharpValence) {
            case 0:
            case 1:
                for (int i = 0; i < N; i++) {
                    faceAverage += faces[i]->centroid();
                    edgeAverage += edges[i]->midpoint();
                }

                faceAverage /= double(N);
                edgeAverage /= double(N);
                nextPos = (faceAverage + 2.0 * edgeAverage + double(N - 3) * pos) / double(N);
                break;
            case 2:
                for (int i = 0; i < N; i++) {
                    if (edges[i]->sharpness) {
                        edgeAverage += edges[i]->midpoint();
                    }
                }

                edgeAverage /= sharpValence;
                nextPos = 0.5 * (edgeAverage + pos);
                break;
            default:
                nextPos = pos;
                break;
        }

        ccVertex = v;
        ccVertex->pos = nextPos;
    }

    if (sharpness > 0) {
        ccVertex->sharpness = sharpness - 1;
    } else {
        ccVertex->sharpness = sharpness;
    }

    ccVertex->opt = opt;
}

void Vertex::setCCLimitVertex(Vertex* v) {
    Eigenstructure eigen = eigenstructures[N-3];

    double pointWeight = eigen.getVectorInverse(0, 0);
    double adjacentWeight = eigen.getVectorInverse(0, 1);
    double diagonalWeight = eigen.getVectorInverse(0, 2);

    Vector3d nextPos = pointWeight * pos;

    for (int i = 0; i < N; i++) {
        nextPos += adjacentWeight * edges[i]->getOpposite(this)->pos;
        nextPos += diagonalWeight * faces[i]->getDiagonal(this)->pos;
    }

    ccVertex = v;
    ccVertex->pos = nextPos;
    ccVertex->sharpness = 0;
    ccVertex->opt = opt;
}

Edge* Vertex::getEdge(int a, int b) const {
    for (int i = 0; i < N; i++) {
        if (edges[i]->equalEndpoints(a, b)) return edges[i];
    }
    return 0;
}

void Vertex::addEdge(Edge* edge) {
    N++;
    edges.push_back(edge);
}

void Vertex::addFace(Face* face) {
    faces.push_back(face);
}

Vector3d Vertex::project(int iK) const {
    Eigenstructure eigen = eigenstructures[N-3];
    Vector3d result = eigen.getVectorInverse(iK, 0) * pos;
    for (int i = 0; i < N; i++) {
        result += eigen.getVectorInverse(iK, 2*i+1) * edges[i]->getOpposite(this)->pos;
        result += eigen.getVectorInverse(iK, 2*i+2) * faces[i]->getDiagonal(this)->pos;
    }
    return result;
}

Vector3d Vertex::derivative(int m, int n) const {
    Eigenstructure eigen = eigenstructures[N-3];

    //iK = 0 is the limit point; it cannot contribute to derivative
    //thus iK = 1, 2 dominate; they have the highest eigenvalues
    Vector3d proj1 = project(1);
    Vector3d proj2 = project(2);
    
    double weight1 = eigen.getx(1, 1, m, n);
    double weight2 = eigen.getx(1, 2, m, n);

    Vector3d result = proj1 * weight1 + proj2 * weight2;

    return result;
}

double Vertex::getCurvature(int start) const {
    if (N != 4) return 1000.0;

    double edgeWeights1[4] = {
         12.0 / 36.0,
          0.0 / 36.0,
        -12.0 / 36.0,
          0.0 / 36.0
    };

    double faceWeights1[4] = {
          3.0 / 36.0,
         -3.0 / 36.0,
         -3.0 / 36.0,
          3.0 / 36.0
    };

    double edgeWeights2[4] = {
         12.0 / 18.0,
         -6.0 / 18.0,
         12.0 / 18.0,
         -6.0 / 18.0
    };

    double faceWeights2[4] = {
          3.0 / 18.0,
          3.0 / 18.0,
          3.0 / 18.0,
          3.0 / 18.0
    };

    Vector3d d1 = Vector3d(0, 0, 0);
    Vector3d d2 = pos * -24.0 / 18.0;

    for (int i = 0; i < 4; i++) {
        int index = (i + start) % 4;
        d1 += edgeWeights1[i] * edges[index]->getOpposite(this)->pos;
        d2 += edgeWeights2[i] * edges[index]->getOpposite(this)->pos;
        d1 += faceWeights1[i] * faces[index]->getDiagonal(this)->pos;
        d2 += faceWeights2[i] * faces[index]->getDiagonal(this)->pos;
    }

    double numer = (d1 ^ d2).magnitude();
    double denom = d1.magnitude();
    denom = denom * denom * denom;

    return numer / denom;
}

void Vertex::calculateCurvatures() {
    double uCurvature = getCurvature(0);
    double vCurvature = getCurvature(1);
    gaussianCurvature = uCurvature * vCurvature;
    meanCurvature = 0.5 * (uCurvature + vCurvature);

    uTangent = derivative(0, 1);
    vTangent = derivative(1, 0);
    normal = uTangent ^ vTangent;
    uTangent.normalize();
    vTangent.normalize();
    normal.normalize();
}

void Vertex::vertexGeometry(bool useMeanCurvature) const {
    double value;
    if (useMeanCurvature) value = meanCurvature * 0.25 + 0.5;
    else value = gaussianCurvature * 0.25 + 0.5;
    glColor3d(value, value, value);
    glVertex(pos);
}

void Vertex::normalGeometry(double length) const {
    glColor3f(0.0, 1.0, 0.0);
    glVertex(pos);
    glVertex(pos + normal * length);

    glColor3f(0.0, 0.0, 1.0);
    glVertex(pos);
    glVertex(pos + uTangent * length);

    glColor3f(1.0, 0.0, 0.0);
    glVertex(pos);
    glVertex(pos + vTangent * length);
}

double Vertex::approxArea() const {
    double result = 0.0;
    for (int i = 0; i < N * 2; i++) {
        const Vector3d& v0 = this->pos;
        const Vector3d& v1 = edges[((i+1)/2) % N]->getOpposite(this)->pos;
        const Vector3d& v2 = faces[i/2]->getDiagonal(this)->pos;
        double a = (v0 - v1).magnitude();
        double b = (v0 - v2).magnitude();
        double c = (v1 - v2).magnitude();
        double s = (a + b + c) * 0.5;
        result += sqrt(s * (s - a) * (s - b) * (s - c));
    }
    return result;
}

double Vertex::approxLocalMeanCurvatureSq() {
    if (N != 4) return 1000000.0;

    double result, temp;

    temp = getMeanCurvature();
    result = 16.0 * temp * temp;
    
    for (int i = 0; i < 4; i++) {
        temp = edges[i]->getOpposite(this)->getMeanCurvature();
        result += 4.0 * temp * temp;
        temp = faces[i]->getDiagonal(this)->getMeanCurvature();
        result += 1.0 * temp * temp;
    }

    result *= approxArea();
    return result;
}

double Vertex::getMeanCurvature() {
    if (N != 4) return 1000.0;
    double uCurvature = getCurvature(0);
    double vCurvature = getCurvature(1);
    meanCurvature = 0.5 * (uCurvature + vCurvature);
    return meanCurvature;
}